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A double well loaded with bosonic atoms represents an ideal candidate to simulate some of the most inter¬ 
esting aspects in the phenomenology of thermalisation and equilibration. Here we report an exhaustive analysis 
of the dynamics and steady state properties of such a system locally in contact with different temperature reser¬ 
voirs. We show that thermalisation only occurs ‘accidentally’. We further examine the nonclassical features and 
energy fluxes implied by the dynamics of the double-well system, thus exploring its finite-time thermodynamics 
in relation to the settlement of nonclassical correlations between the wells. 


The high degree of control available when dealing with ul¬ 
tracold atomic samples makes them ideal candidates for re¬ 
alising prototypical quantum technology devices [1, 2]. The 
range of practical applications that can be addressed using 
platforms based on the physics of ultracold atomic ensem¬ 
bles ranges from metrology and sensing to the achievement 
of quantum memories [3], from ultra-stable atomic clocks [4] 
to the simulation of difficult condensed-matter physics prob¬ 
lems [5], Recently, such range has been extended to quantum 
thermometry [6], while theoretical and experimental interest 
is emerging in the design and implementation of thermody¬ 
namic processes and (elementary) engines based on such sys¬ 
tems [7]. The tuneable interactions among the elementary 
constituents of a cold-atom system, and the availability of 
effective ways of arranging non-equilibrium states of atomic 
systems confined in external optical potentials, provide an al¬ 
most ideal scenario for the study and harnessing of thermody¬ 
namically relevant questions and tasks, indeed recently ther¬ 
mal and number fluctuations have been studied for ultracold 
atoms in two mode traps [8]. 

For such endeavours to succeed, it is absolutely crucial to 
identify a suitable configuration to act as the basic building 
block for a thermodynamic device, and characterise its work¬ 
ing principles in terms of fundamental quantities (such as heat 
and work), which will pave the way to the actual construction 
of the machine itself. 

In this paper, we move exactly along these lines: Inspired 
by the experimental set up of Refs. [7], where a cold atomic 
system is placed in contact with two different thermal reser¬ 
voirs, we consider a slight modification in which the gate po¬ 
tential separating the two reservoirs is replaced by a double 
well potential loaded with a Bose-Einstein condensate (BEC), 
itself a system of vast experimental implementability [9, 10]. 
We set and study explicitly its non-equilibrium dynamics. 
By assuming each well is initially thermalised to its own lo¬ 
cal reservoir, we will show that, in the tunnelling dominated 
regime, a temperature imbalance between the wells leads to 
the emergence of non-classicality, and study how this is linked 
to the equilibration dynamics of the atomic system. Remark¬ 
ably, we show that the genuinely quantum nature of the state 
of the double well does not appear to affect the rate of equili¬ 
bration of the open system at hand. By working in the weak 
coupling regime between each well and its reservoir, which 
allows us to identify clearly the contributions of each well to 
the total heat flux into/out of the local environments, we high¬ 


light a rather rich and complex dynamics of the heat exchange 
across the wells. Further, we examine its relation with the 
emergence of nonclassical correlations within the state of the 
atomic ensemble within a vast range of operating conditions. 


I. DESCRIPTION OF THE MODEL 


We are interested in studying the out of equilibrium dy¬ 
namics and steady-state properties of a system of cold atoms 
loaded in a double-well potential and subject to the effects of 
two reservoirs at different energy. Our Hamiltonian is the two- 
site Bose-Hubbard model [9] given by 77 = 77/ + 77), + 77, 
with [we assume units such that fi — 1 throughout] 


77/ - u >! + -j + a>2 ^ ) > 

OJ U M2-2 , -t2-2\ 

77.9, = — [a\ a l + a' 2 a 2 ), 

77, = -J (a\a2 + aiaty. 


(1) 


Here 77/ describes the free evolution of the atomic systems in 
the two wells, each occurring at the rate set by the single¬ 
atom energy a)j, and with dj, cT the associated annihila¬ 
tion and creation operators for each well. The Hamiltonian 
term 77 s , accounts for the self-interaction (at rate U) between 
atoms occupying the same well, while 77, stands for the tun¬ 
nelling term, which occurs at rate fj. We will focus mostly 
on the tunnelling-dominated regime associated with U - 0. 
However, the interaction-dominated regime corresponding to 
ST — 0, and the intermediate regime will also be addressed. 
The focus of our investigation will be the phenomenology of 
thermalisation of the system, both at the single and two-well 
level. We remark that model, Eq. (1), can be realised in a 
variety of settings including superconducting Josephson junc¬ 
tions [11], trapped ions [12], bimodal optical cavities [13] and 
optomechanical setups [14]. 

While important insight will be gathered by addressing the 
unitary evolution induced by considering 77, the overarching 
goal of this work is the study of the open-system evolution 
created by the contact of the two wells with their respective 
reservoirs. We are interested in addressing the dynamics in- 
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duced by the master equation [15] 

2 

& = -i[#,e,] + 2^te»). (2) 

7=1 

where we have introduced the overall-system density matrix 
at a generic time t, g r , and the Lindblad super-operators 

-Cj(Qt) =7j(nj + \ )ia j g t a j - ^{ata,-, £>,}) + 

' ' (3) 

yjrijU\g t aj - 

which describe the incoherent particle-exchange process 
(occurring at rate yj) between a well and the respective 
reservoir (assumed to have a thermal occupation number nj). 
Eq. (2) is the key equation in our analysis to follow. We 
remark that in certain working conditions this description 
of the dynamics is not always valid. In particular, when 
the scattering length of the BEC is large, non-Markovian 
dynamics can play an important role [16]. We therefore 
assume that the scattering length is sufficiently small to 
ensure the validity of the Markovian approximation [16]. 


II. EXACT SOLUTIONS OF THE 
TUNNELING-DOMINATED REGIME. 


mixing transformation Uj(9 ) = exp[-i|(aja 2 + «i a^)] with 
9 = -\ arctan (27/A), which leaves us with the new quadratic 
Hamiltonian 

( H q /io i = Oi(X[ + p\) + fl 2 (xi + pi), (7) 

describing two freely evolving harmonic oscillators at the re¬ 
spective frequencies fli = 1 + (A - E)/2, Q 2 = 1 + (A + E)/2, 
with T = VA 2 + 47 2 . For a Gaussian initial state of the sys¬ 
tem [18], given the quadratic nature of the Hamiltonian, rather 
than tracking the evolution of the density matrix of the system, 
we can restrict our attention to the evolved form of the covari¬ 
ance matrix <r of entries cr,y = ({Pi, Pj}) - (Pj)(Pj), where 
Pi’s are the elements of the vector of quadrature operators 
P T = (x\ pi x 2 pi) and the expectation value of such vec¬ 
tor (calculated over the state of the system), which bear full 
information on the state of the system. Both are readily gath¬ 
ered as 

cr u (t) = Mcr „(0)M T , (P) t = M(P) 0 , (8) 

with M = T(9)R\(D.it)R 2 (Q. 2 t)T(0) T , cr„(0) [(/%] the covari¬ 
ance matrix [vector of phase-space displacements] of the ini¬ 
tial state of the system and cr u (t) [(H)/] its time-evolved ver¬ 
sion. In Eq. (8) Rj(£ljt) (j =1,2) and T(9) are the symplectic 

transformations corresponding to the free evolution e ' n J a j a J t 
and two-mode mixing Ut(9). Explicitly 


In order to gather insight into the basic coherent processes 
of the system in the case of tunnelling-dominated regimes, 
we set U - 0 in 'PI and address the unitary evolution first. 
We define the canonical quadrature operators {xi, p\, x 2 , p 2 ) 
as [17, 18] 

*; = ^( a ; + a D’ and = (4 > 

and recast the Hamiltonian into the form 

= Y (x 2 i + p\ + l) + Y (*2 + Pi + l) -J (xix 2 + Plpl) , 

(5) 

with 1 the identity operator. By neglecting trivial constant 
terms, Eq. (5) can be thus interpreted as a quadratic form iden¬ 
tified by the adjacency matrix 


(CD 1 0 -J 0 ) 

0 a>i 0 -J 
-J 0 cu 2 0 
, 0 ~sr 0 u> 2 , 


( 6 ) 


which has been written in the ordered operator basis 
{x\,pi,X 2 ,p 2 ). In what follows, we rescale all the relevant 
frequencies with respect to a>\. In these units, we have 
W 2 —» W 2 /«i = 1 + A with A a dimensionless bias between the 
two wells, and 77 —> 7 = J'/oj]. The rescaled Hamiltonian 
77/w 1 can be diagonalised by means of a simple two-mode 


Rj(Cljt) 


I cosfO jt) sin(G 7 f)\ 
sin(D jt) cos(QjO ) ’ 


cos 9 

0 

sin# 

0 

0 

cos 9 

0 

sin# 

- sin# 

0 

cos # 

0 

0 

- sin# 

0 

cos # 


(9) 


We now concentrate on the situation where the particles in 
each well are initially at thermal equilibrium with their local 
reservoirs. The initial covariance matrix will thus be that of a 
two-mode thermal state 


' 2«i +1 0 0 O' 

0 2Ti\ + 1 0 0 

0 0 2n 2 + 1 0 

„ 0 0 0 277 2 + 1 j 


( 10 ) 


with rij = (a/ij) the mean number of particles in the / lh well. 
For such an initial state, the phase-space displacements are all 
null and full information on the evolved state is provided by 
the covariance matrix 


o- u (t) = 


n 1 
0 
Cl 

c 2 


0 c\ ci’ 

III -Cl Cl 

-ci n 2 0 

ci 0 n 2 , 


( 11 ) 
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with elements 


4/A(«i - ht) sin (Tf/2) 2J (n\ -Hi)sin(Tf) 
ci = -^-, c 2 = - 


«i = 


T 1 r 

4/ 2 («i - h 2 ) cos (rr) + 4/ 2 (77 1 + h 2 + 1) + A 2 (2h! + 1) 


4/ 2 (hi - n \)cos (r r) + 4J 2 (ni + h 2 + 1) + A 2 (2hi + 1) 

n 2 = --- f2 - 1 - 1 -• 

( 12 ) 

If both wells are at the same initial temperature, i.e. n\ = n 2 , 
then c\ = C 2 = 0 and n\ = n 2 = 2n\ + 1, i.e. the system does 
not evolve in time and the two wells remain at their thermal 
equilibrium, notwithstanding the tunneling. This is a clear 
interference effect. Moreover, for identical single-atom en¬ 
ergy in each well, i.e. A = 0, ci is null, showing that the 
position [momentum] xi [p\] gets correlated with p 2 [x 2 ], In 
general, such correlations do not imply necessarily the setting 
of entanglement between the wells [19], Indeed, the tunnel¬ 
ing term of the Hamiltonian, < H, in Eq. (1), can generate en¬ 
tanglement only when the state of at least one of the wells 
is sufficiently non-classical. In the context of our investiga¬ 
tion here, this basically implies the preparation of squeezed 
states of the wells [20]. This can be understood by noticing 
the formal analogy between 7Y, and the generator of a two¬ 
mode mixing transformation and considering, for the sake of 
argument, the resonant case A = 0. Under such conditions, 
moving to the interaction picture with respect to *Hf, the time 
evolution operator would correspond to Uj(Jt), which gives 
rise to no entanglement between the two wells when they are 
prepared in thermal states (even at different effective temper¬ 


atures), as demonstrated in Ref. [20], However, this does not 
imply that the dynamics of the two-well system is trivial. In 
fact, in general, quantum correlations (of a form weaker than 
entanglement) are set by U(Jt) when acting on thermal states 
with «i 4 h 2 . We will address the emergence of discord-like 
quantum correlations [21-23] and its relation to the inter-well 
exchange process in a later section. 

We now move to solving the full dissipative dynamics gov¬ 
erned by Eq. (2) for U — 0. The problem can be effi¬ 
ciently solved by using a suitable Gaussian ansatz: We first 
translate Eq. (2) into the phase space by deriving a differ¬ 
ential equation for the symmetrically ordered characteristic 
function x(Pi,P 2 , t) = Ti[Di(J3i) ® D 2 iP 2 )Qt\ [17, 18]. Here 
Dj(/3j) = exp \J3ja\ -P*cij] is the Weyl displacement operator 
with amplitude f5j 6 C for system j = 1,2. Using the phase- 
space relations [18] 


^Dj(Pj) ~ + ^Dj(J3j), 

IB a (13) 

DjiPjTa ~ (~2 - 

after a lengthy but otherwise straightforward calculation, we 
find that Eq. (2) takes the form of the Fokker-Planck equation 


d t X<J3x 




_d_ 

Wi 



_d_ 

+ /3l d/?2 


i 

j= i 


CJi 


d_ _ d_ 

P j dPj P'W 


j J 


Zl 
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By letting y3 ; = xj + ipj [so that x(J3uPi,t) -» 

X(xi,pi,x 2 ,p 2 ,t)] and expressing the characteristic function 
in terms of the entries of the vector of quadrature variables, 
we can write x( x i,pi,x 2 ,p 2 ,t) = exp[/P T A - \P T &P], where 
we have introduced the generic vector of C-numbers X T = 
(yi z,\ y 2 z 2 ) and matrix & whose elements we aim at finding, 
which we do by solving Eq. (14). In Methods we provide 
the set of differential equations for the elements of X and <x 
obtained when evaluating both sides of Eq. (14) and equating 
them term by term. 

The explicit solution of the problem at hand leads to a time- 
evolved covariance matrix of the general block form 


cr(0 = 


(m\ 1 

U T 



(15) 



Pi 

\ 


_d_ 

Wj 



+ Tj 




XiPuPi)- 


(14) 


where c is a 2 x 2 matrix of correlations among the quadra¬ 
ture operators of the system. The diagonal structure of the 
blocks pertaining to the the individual wells shows that, lo¬ 
cally, the system thermalises at temperatures determined by 
the explicit form of m\ 2 - However, as c is, in general, not null, 
global thermalisation is not achieved: the overall system never 
thermalises, notwithstanding an explicitly dissipative evolu¬ 
tion. This is clearly seen by looking at the general form of 
the steady state. Although the analytic form of the non-zero 
elements is readily achievable for any value of the parame¬ 
ters involved in the problem, they are, in general, too cumber¬ 
some to be reported here. However, assuming y i = y 2 = y. 
the steady-state of the system is determined by the covariance 
matrix 
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0~ss — 


47 2 (/ii+«2+l) 

y 2 + A 2 


(2ni + 1 ) 0 

q 47 2 (ni +n 2 + 1) 

2/A(wi~W2) 


y 2 + A 2 
2Jy(Jj\-n 2 ) 
y 2 + A 2 


+ ( 2 n, + 1 ) 

27y(ni-w 2 ) 

■y 2 +A 2 
27A(«j -«2) 
y 2 +A 2 


_ 2JA(n\-n 2 ) 
y 2 + A 2 
2Jy(n\-n 2 ) 
y 2 + A 2 

4 j 2 p'J + 1 ) +( 2 « 2 - 

0 


2Jy{n\-n 2 ) 
y 2 + A 2 

_ 2JA(7?i -m) 
y 2 +A 2 


■i) 


o 


4/ 2 (n[+W2+l) 


- 4 - 


(16) 


(a) (b) 



FIG. 1: (a) Maximum fidelity between the instantaneous state of the 
system and a globally thermal state, plotted against the dimensionless 
evolution time and the (dimensionless) bias between the energies of 
the wells A. (b) Corresponding estimate of the mean energy p of the 
target globally thermal state. In both panels, we have taken 7 = 2, 
= zF /2 = 1 . 

with £ = 4 J 2 +y2+A 2 • Clearly, only for n\ = n 2 the structure 
of the global covariance matrix takes a thermal-like form. 
However, this does not preclude the possibility to achieve 
accidental thermalisation, i.e. situations such that the state of 
the system either becomes globally/locally thermal, or closely 
approximates an equilibrium configuration. This will be the 
focus of the following analysis. 


III. ASSESSMENT OF DYNAMICAL THERMALISATION 


We shall start with the study of the unitary case. Dynam¬ 
ical thermalisation in closed-system dynamics is a topic of 
vast interest, which has recently attracted considerable atten¬ 
tion at both the theoretical and experimental level [24]. Our 
approach is based on the assessment of the distance between 
the time-dependent state of the system and a generic (either 
global or local) thermal state. Quantitatively, as a measure 
of the distance between two states p\ 2 , we use the Ulhmann 
fidelity [25] 


F(pi,p 2 ) = ^Tr ^yfplpi y[p~\} ■ (17) 

For Gaussian states, it can be conveniently evaluated using 
the covariance matrices <-r i2 associated to the states under 
scrutiny. The explicit formula, which has been recently re¬ 
ported in Ref. [26], reads 


Fieri, o- 2 ) = 4 


(yfx + Va' - l ) 2 
Vdet(cr! + oq) 


(18) 


where Q = z(cr v © cr y ) is the two-mode symplectic matrix ( cr y 
being the y-Pauli matrix) and x — 2 V T~\ + 2 VJ 2 + 1/2 with 

^ detfGoqGcri - I 4 ) det(<xi + zn)det(cr 2 + z'Q) 

16det(<xi + 0 * 2 ) ’ ” 16det(<n + cr 2 ) 

(19) 

In our case we consider <r\ = cr u (t) [cf. Eq. (11)] and <r 2 
given by either cr^ = ( 2 p + OI 4 , i.e. the covariance matrix of 
a globally thermal state with mean number of excitations p, or 
cr/- = ( 2 pi + 1 )H 2 © ( 2 p 2 + 1 ) 1 - 2 , which is the one associated 
with the tensor product of locally thermal states (each with 
mean number of excitations pj). For clarity, we have indicated 
with II ( , the identity matrix of dimension n. 

We present the case of global thermalisation first: after 
calculating the time behaviour of F(cr u (t ), cF,') for various 
choices of A, we have numerically evaluated the value of p 
that achieves the maximum of F(cr u (t), cr ^). In Fig. 1 we 
show both such value and the corresponding estimate for p. 
The state fidelity remains evidently quite large, being only 
partially depleted by an increasing value of A (the depen¬ 
dence on such parameter is quite non-trivial, given that for 
A = 2.5, for instance, values very close to those associated 
with A = 0 can be achieved, at suitable times in the evolution). 
However, while at small values of A the target state changes 
very little with time, this is not the case for increasing bias: 
the value of p corresponding to a non-zero A oscillates with 
a non-negligible amplitude as this parameter grows. In any 
case, perfect thermalisation is never achieved, a result that is 
strengthened by the analysis that we will report in the next 
Section. 

The situation is somehow different when locally thermal 
target states are considered [cf. Fig. 2]: besides the expected 
times at which a full period of the evolution is achieved, it is 
possible to identify instants of time at which the state of the 
double-well system is indeed very close to a locally thermal 
state ( F(cr u (t ), cri;) > 0.999), which would suggest the occur¬ 
rence of accidental dynamical thermalisation. 

In the open-system dynamics case, a similar calculation al¬ 
lows us to evaluate the fidelity with both a globally thermal 
and locally thermal state as shown in Fig. 3, which studies 
the effects of both the energy bias [panel (a)] and a differ¬ 
ence in the damping rates of the two wells [panel (b)]. Quite 
evidently, both effects spoil the state fidelity, which however 
achieves values that are either precisely 1 or very close to it. 
Indeed, focusing on the unbiased case with A = 0 and y 1 = y 2 , 
we know that the solution is given in the form of Eq. (15). 
Moreover, the off-diagonal block matrix c turns out to be anti¬ 
diagonal with entries equal in modulus but opposite in sign. 
Therefore, in order to determine if the system has accidentally 
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(a) 



FIG. 2: (a) Maximum fidelity between the instantaneous state of the 
system and a locally thermal state, plotted against the dimensionless 
evolution time and the (dimensionless) bias between the energies of 
the wells A. (b) & (c) Estimate of the corresponding mean number 
of excitations 2 of the target locally thermal state. In both panels, 
we have taken J -2, n\ = t? 2 /2 = 1. 

thermalised, we need only to determine if, at some value of t, 
these entries are identically zero. After some manipulation. 


(a) 

F(a(t),o G 2 ) 






FIG. 3: Open-system dynamics, (a) Fidelity with a globally thermal 
state for J = 2, ri\ = 1, n 2 = 2, yi/tui = y 2 /w, = 1- We have taken 
A = 0 (blue line) and A = 0.5 (red line), (b) Fidelity with locally 
thermal states for J = 2, = 1, n 2 = 2, and A = 0. We have taken 

7 i/o)] = y 2 /wj = 1 (red line) and yi/wi = 1 with y 2 /o)j = 3 (blue 
line). 


this condition reduces to the transcendental equation 

2T 

e ri ' = cos (2 Jt) - — sin (2 Jt) . (20) 

n 

Interestingly, this ‘accidental’ thermalisation is independent 
of the temperature of either well and only concerned with 
the tunnelling strength and the damping rate. For the same 
parameters taken to obtain the red curve in Fig. 3 (b), 
we find Eq. (20) has two solutions: t ~ 1 .03438oj 1 ! and 
t ~ 1.33749oq 1 , clearly corresponding to the two instances 
of local thermalisation in Fig. 3 (b). Furthermore, we find the 
thermal occupation numbers of the wells at the first instance 
of thermalisation are Ti\ = 1.597 and n 2 = 1.403, and at the 
second are Ti\ = 1.422 and «2 = 1.578, thus suggesting that 
the two instants of accidental thermalisation correspond to 
an almost swap of the two local thermal states. Increasing J 
leads to more instances of accidental thermalisation occurring 
before the system equilibrates to its steady state. 


IV. ASSESSMENT OF THE NON-CLASSICAL NATURE OF 
THE STATE OF THE SYSTEM. 

Values of fidelity so close to unity should not lead to mis¬ 
interpretation of the actual nature of the state of the two-well 
system. In fact, any assessment of fidelity should be accom¬ 
panied by the study of problem-specific figures of merit able 
to provide a more fine-grained characterisation of the state at 
hand. For the sake of a study on thermalisation, a significant 
class of such quantifiers is embodied by measures of quantum 
correlations. 

In this respect, it is important here to assess the role, if 
any, various forms of quantum correlations play in the dy¬ 
namics highlighted above. It is quickly confirmed that, as an¬ 
ticipated before, the system never becomes entangled. While 
this is expected in light of the nature of the interaction and ini¬ 
tial state being considered, nothing prevents the settlement of 
weaker forms of quantum correlations, such as quantum dis¬ 
cord (QD) [23]. QD is the difference between two classically 
equivalent definitions of mutual information when applied to 
a quantum system [21, 22]. A non-zero degree of QD implies 
that, in a bipartite system composed of parties A and B, infor¬ 
mation can be gathered on system A by interrogating party B. 
For Gaussian states, QD is captured by the Gaussian quantum 
discord [27-29], which entails that the interrogation of B only 
involves Gaussian measurements. For a generic covariance 

matrix S = | ^ j, QD is then defined following Ref. [27] 

(to be consistent with the definition of the vacuum state used 
throughout) 


with 


D c =h(y[T l )-h(d-)-h(d + ) + h(VE^), (21) 













6 


E mm = 


1 


(A - l) 2 

1 


2l\ + 


(/i - l)(/ 4 - A) + 2|/ 3 | ^/f 2 + (A - l)(/ 4 - A) 


for (7 4 - A A) ^ A (A + A)(A + l)' 


2 /, 


A A - 7 2 + / 4 - ^ + (/ 4 - /i/ 2 ) 2 - 2/|(/ 1 / 2 + / 4 ) 


otherwise, 


where 


present. Also we see the system tends to equilibrate faster. 


h(x) = 



dl= l -(A± VA 2 -4A), 


A = A + A + 2/ 3 , 



( 22 ) 


and 7] = detA, I 2 = detB, 7 3 = detC, and 7 4 = detS. In 
Fig. 4 (a) we study QD against the energy bias for the case of 
the unitary solution Eq. (11). Intuitively we would expect that 
for A = 0, owing to the full symmetry enforced in the system, 
QD will be maximised. This is indeed the case, as it can be 
seen in Fig. 4 (a). However, an interesting feature appears as 
we increase the bias. At A/J ~ 2.5, n\ = 1, and Th = 5, QD 
exhibits a plateau, which implies the existence of an ‘optimal’ 
value of the bias, dependent on the temperature difference, 
that helps amplify the non-classicality of the system. Further 
increase of A pushes the systems too far off resonance, and the 
coherence decays. In Fig. 4 (b) we examine this phenomenon 
closer, for a fixed temperature difference and small biasing, 
A/J — 1, (red), we see the oscillatory behaviour changes and 
the first zero-point is lifted. At the optimal value of A (solid 
black) the plateau is clearly evident. When we increase the 
bias further, we see the decay in the non-classicality, as well 
as a change in the periodicity of the system. 

Turning our attention to the dissipative case, in Fig. 5 we 
compare the unitary dynamics with the dissipative case for un¬ 
biased wells [panel (a)], and biased ones [panel (b)] at various 
differences in temperature. For unbiased wells, we see dissi¬ 
pation quickly suppresses the the oscillations and we reach a 
steady state with non-vanishing QD. As we increase the tem¬ 
perature between the wells we see an increase in the QD for 
the unitary case and the steady state QD is larger for increas¬ 
ing temperature difference. When we bias the wells, taking 
A/J = 2.5 for all temperature differences, we see the dissi¬ 
pative dynamics clearly show the enhanced non-classicality. 
While the time to reach equilibrium appears unaffected, the 
steady state is significantly more non-classical than in the un¬ 
biased situation. This may imply that in this situation the non- 
classicality plays no role in reaching equilibrium. In Fig. 5 
(c) we examine the effect that self-interaction has on the dy¬ 
namics of nonclassicality. In order to do so, we compute the 
Gaussian discord of the hypothetical Gaussian state having, as 
covariance matrix, the one achieved by calculating the entries 
cTij over the non-Gaussian state resulting from a chosen non¬ 
zero values of U. Evidently the larger the self interaction, the 
more self ordered each well becomes, diminishing the effect 
of the tunnelling and reducing the amount of nonclassicality 


The nonclassicality of the steady state is delicately depen¬ 
dent on the temperature difference, as well as the tunneling 
rate and the bias. In Fig. 6 we examine this behaviour closer, 
fixing Ti\ and y\/u)\ = 72/^1 = 1 with J = 2 and A = 5. 
The only conditions for which the system does not exhibit 
nonclassical correlations is the trivial one of = «i • As we 
increase the temperature imbalance we see that QD increases 
to a maximum value before slowly decaying [cf. Figs. 6 (a) 
and (e)]. If we fix the temperature difference such that 77 1 = 1 
and «2 = 5, we see in panels (b) and (c) that there are optimal 
values of the remaining parameters that give the largest value 
of discord. While the reservoirs have been kept at moderately 
low energies, in panel (f) we significantly increase both 77 1 
and 77 2 and see that large values of QD can still be achieved. 
An unbiased configuration leads to values of discord of the 
order of 10 \ Increasing the bias, such values are raised by 
up to one order of magnitude. 


(a) 



(b) 

£>G 



FIG. 4: (a) Discord versus bias and evolution time for the case of 
closed-system dynamics. We have taken = 1, n 2 = 5, J = 2. 
(b) Behaviour of quantum discord in the open-system scenario for 
= 1, n 2 = 5, J = 2 and the bias choices A = 1 (red curve), 5 
(black curve), and 10 (blue curve). 
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(a) 

D c 





FIG. 5: (a) Dynamical discord for the unitary (gray) and dissipative 
(red) cases. For both 7 = 2, H, = 1, and «2 = 2 (dotted), 5 (dashed), 
and 10 (solid). For all the dissipative cases y\!u>\ = y 2 /(±>i = 1. 
(b) As for panel (a) but with the optimal bias (for 77i = 1, no = 5) 
between the wells, A = 5. (c) Dissipative dynamical discord 7 = 2, 
A = 4, yi/wi = 72 /wi = 1, Hi = 1, and n 2 = 2, with self interaction 
term U = 0, U = 1 and (7 = 3 going from top to bottom. 


V. DYNAMICS OF THE ENERGY FLUX BETWEEN THE 
WELLS 


It is important to gather insight into the details of the ex¬ 
change of energy between the wells of the system, which is 
at the basis of the process of quasi-thermalisation highlighted 
so far and takes place in two forms: an exchange of particles 
between the wells and a similar process occurring at the inter¬ 
face between the double-well system and the reservoirs. The 
aim of this section is to identify the contribution coming from 
both such fluxes. We are thus interested in quantifying the flux 
into/from well j — 1 , 2 , which we label as < 2 , , and the total 
flux <2, of . These are given by the quantities 

<2,„, = Tr[77<3,p], Qj = Tr['H J d,g j ], (j = 1,2), (23) 

where 77 ; = ojjia .'dj + 1 /2) is the free evolution of a single 
well and Qj is the density matrix of well j. Conveniently, these 
quantities can be directly evaluated from the covariance ma¬ 
trix (and we will assume both wells to have the same damping 



FIG. 6 : Steady-state discord between the two wells for yi/a>i = 
72 /wi = 1, and Hi = 1. (a) Plotted against H 2 for 7 = 2, A = 5. 
(b) Against A for 7 = 2 and H 2 = 5. (c) Against 7 for H 2 = 5 
and A = 5. (d) Maximum discord attainable for a given value of n 2 
found by optimising with respect to both A and 7 when Hi = 1 and 
yi/rn, = y 2 /mi = 1. (e) Steady-state discord studied against both Hi 
and H 2 when 7 = 2, A = 5, and 71 /tui = y 2 /wi = 1. The black line 
shows that Do is identically null only when Hi = n 2 . (f) Steady-state 
discord against n 2 for Hi = 100 with 7 = 2, yi/w, = y 2 /a>i = 1. 


rate, i.e. y\/u>\ = 72/01 = 7 ). We find 


2/ 2 (l + A)(H 2 - HO . 

sin 


Qi = e-* ._ 

V47 2 + A 2 

<2 2 = -(l + A)<2i, <2, 0 , = 0 . 


( V47 2 + A 2 t) 


(24) 


From here it is easy to confirm that <2 1 /w 1 = -<2 2 /o 2 and 
clearly taking 7 = 0 or A = 0 recovers the unitary and unbi¬ 
ased limits respectively. However, this behaviour is only for 
the special case of the wells initially being thermalised with 
their baths, while taking a different initial state this behaviour 
no longer holds. Indeed, what is special about our initial state 
is that it conserves the total energy of the system. This is read¬ 
ily seen given that < 2 , 0 , = 0 and it is easy to confirm that 

(2 tot = 1 + ^ + «i + (1 + A)n 2 , (25) 

for all t and J. Of course, the energy of the individual 
wells changes dynamically (until settling into the same steady 
state). 

We can gain further insight into the reason for this by ex- 
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amining closer the quantity we are calculating, i.e. 


0,0, = Tr VH3 tQ ] 


= Tr 





= -iYv\ r H( r ftp- ,yfi)\ + '\v 




i=i 



The tunnelling term in Eq. (1) commutes with and when 
U - 0 the only contribution to the total flux is from the free 
evolution of each well. Therefore we are interested in calcu¬ 
lating 

Q, ol = Tr [(hi + n 2 ) ^(g)] + Tr |(h, + n 2 ) Zi(g)}. (26) 


In a tedious but otherwise straightforward calculation, we can 
explicitly evaluate this expression when assuming the special 

initial condition g = e ® e , with Z,j = TrEe - ^^']. We 
find that both the terms entering Eq. (26) are identically zero, 
thus showing that, in the U - 0 case, the net heat flux is null 
due to two special circumstances: on one hand our chosen 
initial state, on the other hand the tunnelling term commutes 
with the super operators. 

In Fig. 7 (a) we show the dynamics of the various energy 
fluxes given in Eq. (23). We notice how the flux into the cooler 
well is proportional to the flux out of the hotter well, which 
results in a null net flux. Needless to say, the single-well fluxes 
only account for the net intake/outtake of particles for one of 
the wells and do not provide information on the actual balance 
between the contribution due to the coupling to the reservoir 
and that due to the coherent inter-well interaction. 

We can study the intermediate dynamical regime where 
the self-interaction is non-zero and comparable with the 
tunnelling by numerically solving Eq. (2) and examining the 
behaviour of the heat fluxes, of which we illustrate some 
examples in Fig. 7 (b) and (c) [we refer to the caption for 
an account of the parameters used in the simulations]. The 
total flux is now non-zero, and the energy is not conserved. 
However, the average occupation number is conserved, i.e. 
(a\ai)/u>\ = -{a\a 2 ) / w 2 , which follows directly from the 
previous arguments. 


VI. DISCUSSION 

The analysis above shows that neither global nor local ther- 
malisation with the reservoirs is achieved. The fidelity be¬ 
tween the density matrices of the time-evolved state and the 
target thermal one (whether globally or locally) connects the 
closeness of the populations of the energy levels of the former 
to the statistics of the latter. However, the interaction between 


the wells establishes strong quantum coherence between the 

(a) 

<3 





FIG. 7: Steady-state discord between the two wells. In all panels 
J = 2, y = 1, «i = 1, «2 = 2, A = 4. Red Line: hotter well (well 
2), gray line cooler well (well 1), and blue the total flux. With (a) 
U = 0. (b) U = 1. (c) (7 = 3. 


particles of the systems, which in turn results in the genera¬ 
tion of a substantive degree of quantum correlations, albeit of 
a nature weaker than entanglement, which prevents the ther¬ 
mal character of the resulting state. 

The analysis reported here also has the merit of providing 
rather deep insight into the phenomenology of quantum 
correlations between the wells. We have qualitatively and 
quantitatively examined the dynamics and steady state of a 
BEC loaded into a double well potential. While the wells 
remain separable at all times, thus sharing no entanglement, 
by exploring the behaviour of the quantum discord we 
find the system to be always non-classical, except under 
trivial, uninteresting conditions. Furthermore, the degree 
of nonclassicality of the system is dependent on the energy 
bias between the two wells. For identical wells, a significant 
amount of QD is possible, provided that a large temperature 
imbalance is established. Such nonclassicality can be greatly 
enhanced by taking a suitable value of tunnelling, which must 
be a function of the given bias. The transfer of heat in the 
system is equally complex. 
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VII. METHODS 

Differential Equations. Here we provide the complete set of differential equations that describe the dissipative dynamics 
considered throughout. 


d,yi = -Jz 2 - j-y i + w izi> 


d,cr\ = yi + 2n 1 y 1 - lJcr{ 2 - no-j ■ 
d,cr p - yi + 2Jiiyi + 2Ja<~ px - yicrj’ - 2o>i(Xj 
5 » cr i P = ^ (°1 2 _ cr ^) - 710*1 - "1 (°1 - <) - 


<9,T2 = —Jzi - yy 2 + w 2 Z2, d,z\ = jyi - fzi - u>iyi. 


■ 2io\cr xp , 

r xp 


<9,cr]\, = -J (cr xp + u xp ) 
dto-'l p 2 = J (erf - tr P 2 ) - ^2 


ri +72_v 


cr 


12 


+ a»i cr^ 2 ' + W 20 "j 2 ’ d f crf 2 = / (i 


xp p 

cr ]2 + 6 l»i cr 


12 


■ <u> 2 cr 


12 ’ 


d t cr x 2 - y 2 + 2n 2 y 2 

3,cr£ = y 2 + 2n 2 y 2 + 2,/crj 2 

5,o-f - 


d t zi = Jyi - yZ2 - to 2 y 2 . 


2J(T PX ■ 

xp 


= J (°"l2 - 


cr 




12 

)- 


y 2 cr^ + 2(jj 2 cr x J’, 
y 2 oi - 2 lj 2 ct x 2 p , 
) - 72 of - to 1 (cr- 2 - o|), 


yi+ra p 
2 u 12 


xp 

■ a»icr 12 ■ 


nx 

to 2 o- 12 , 


A rrP X - T (rr x rr P \ - 21221 rr PX - / , rr x i , , 

O t cr l 2 -J\cr 2 cr l J 2 cr 12 Wi<t 12 + a»2cr 12 - 


(27) 


Discussions on self interaction dominated limit. While the 
main analysis treated the tunnelling dominated regime, the 
opposite extreme is determined setting J = 0 and explor¬ 
ing the situation where self-interaction dominates. In this 
instance, the two wells are completely decoupled from one 
another. We can directly solve Eq. (2) by projecting onto 
the number states |n). Since these states are eigenstates of 
the Hamiltonian for J - 0 the steady-state will be entirely 
independent of U. In fact, regardless of the initial state we 
find the steady state for each well to be pj = -2-e~^ iC ° inj 

with — =—fr, which is the Boltzmann distribution for a 

harmonic oscillator with thermal occupation hj. Clearly then, 
if our initial states are already thermalised with their local 
reservoir we see no dynamics. For any other initial state. 


the two wells thermalise independently to their respective 
reservoir temperatures. 
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